program mkdeco;

{Pittsburgh-Ithaca, summer 99}

{
 INPUT FILES:
 decoration.inp
 symmetry.def
 }

uses common,tools;
Const
   listsize = 200;
Type
   i7	       = array[1..7] of integer;
   v2list      = array[1..t_sym] of v2;
   v5list      = array[1..t_sym] of v5;
   i1list      = array[1..listsize] of integer;
   i5list      = array[1..listsize] of i5;
   noderec5    = array[1..listsize] of i7;
   atomreclist = array[1..listsize] of atomrec;
   v5sym       = array[1..t_sym] of v5;
   v22	       = array[1..2] of v2;
   tsyml       = array[1..t_sym] of integer;
Var
   ofla				      : i3;
   b				      : i5;
   u,v,dv,w,v0			      : v2;
   eofile,figout		      : boolean;
   txt,txt2			      : text;
   code				      : REAL;
   i,j,k,jk,l,m,n,o,iat,jat	      : integer;
   dbg,i0,nstar,dps		      : integer;
   jch				      : integer;
   n_pxyz0			      : integer;
   cmp,por,zz			      : myreal;
   ix10				      : i10;
   x10				      : v10;
   x2				      : i2;
   defdp			      : tsyml;
   defdpz			      : array[-1..1] of i2;
   autocc			      : array[1..t_obj] of integer;
   obj_dpz			      : array[1..t_obj] of tsyml;
   obj_syml			      : array[1..t_obj,1..t_sym] of i3;
   obj_flad			      : array[1..t_obj] of i3;
   syml				      : array[1..t_sym] of i3;
   obj_bas_dp			      : array[1..t_sym,1..t_sym] of tsyml;
   obl				      : array[1..t_sym,1..t_sym,1..t_sym] of i2;
   obj_bas_list			      : array[1..t_sym] of i2;
   obj_bas_mu			      : integer;
   obj_bas			      : array[1..2] of i10;
   objori,or2			      : array[1..t_obj,1..t_sym] of v2;
   obj_break			      : array[1..t_obj] of boolean;
   bobor			      : array[1..t_obj] of integer;
   mu_sb			      : array[1..t_orb] of integer;
   ch_sb			      : array[1..t_orb] of tsyml;
   ptr				      : array[1..t_obj] of integer;
   ptra				      : array[1..t_obj] of v2;
   ptra5			      : array[1..t_obj] of v5;
   ptrz				      : array[1..t_obj] of real;
   ncode,io			      : integer;
   liftperplim,lifteps		      : myreal;
   ctrl				      : integer;
   blift			      : boolean;
   bv1,bv2,bv12			      : boolean;
   papo				      : i3;
   pxyz_rsd			      : array[1..3*t_orb] of myreal; {positional parameters}
   ss				      : mystring;
   ntobj,ntorb,nu_offs,nmir,n2s,nfold : integer;
   first_atom			      : boolean;
   ob				      : tsyml;
   nfli				      : integer;
   fli				      : array[1..t_sym] of i3;
   or5				      : v5sym;
   ocolor			      : integer;

procedure v10v2(b1,b2 : i10;x:v10; var y:v2);
var i,o	: integer;
   bb	: v22;
   v	: v2;
begin
   for o:=1 to 2 do v[o]:=0.0;
   for o:=1 to 2 do
      for i:=1 to 10 do v[o]:=v[o]+x[i]*re[i,o];
   for i:=1 to 10 do if (b1[i]>0) then bb[1]:=re[i];
   for i:=1 to 10 do if (b2[i]>0) then bb[2]:=re[i];   
   y[1]:=(v[1]*bb[2,2]-v[2]*bb[2,1])/(bb[1,1]*bb[2,2]-bb[1,2]*bb[2,1]);
   y[2]:=(v[2]*bb[1,1]-v[1]*bb[1,2])/(bb[1,1]*bb[2,2]-bb[1,2]*bb[2,1]);
end; { v10v2 }
procedure r5r2(b2 : i2;x:v5; var y:v2);
var i,o	: integer;
   bb	: v22;
   v	: v2;
begin
   {NOTE: e[i] is 5-fold set, re[i] is 10-fold set!}
   for o:=1 to 2 do v[o]:=0.0;
   for o:=1 to 2 do
      for i:=1 to 5 do v[o]:=v[o]+x[i]*e[i,o];
   for i:=1 to 2 do bb[i]:=re[b2[i]];
   y[1]:=(v[1]*bb[2,2]-v[2]*bb[2,1])/(bb[1,1]*bb[2,2]-bb[1,2]*bb[2,1]);
   y[2]:=(v[2]*bb[1,1]-v[1]*bb[1,2])/(bb[1,1]*bb[2,2]-bb[1,2]*bb[2,1]);
end; { r5r2 }
procedure orb_union(mu0	:integer; Var mu : integer; Var r:v2list; Var w:v5list; Var x:atomreclist); {in read_raw_deco}
var i,j,o : integer;
   b	  : array[1..t_sym*t_sym] of boolean;
   eps,ar : myreal;
   v	  : v2;
begin
   eps:=0.0001;
   for i:=mu0 to mu do begin
      b[i]:=true;
      {write('IN_orb_uni ',i:1,' ');
      for o:=1 to 2 do write(r[i,o]:1:4,' ');
      writeln;}
   end;
   for i:=mu0 to mu-1 do if b[i] then
      for j:=i+1 to mu do begin if b[j] then
	 for o:=1 to 2 do v[o]:=r[i,o]-r[j,o];
	 ar:=sqrt(v[1]*v[1]+v[2]*v[2]);
	 if (ar<eps) then b[j]:=false;
      end;
   j:=mu0-1; for i:=mu0 to mu do if b[i] then begin
      j:=j+1; r[j]:=r[i]; w[j]:=w[i]; x[j]:=x[i];
   end;
   mu:=j;
end; { orb_union }

procedure obj_union2(no,center : integer;Var gor:v5sym); {in read_raw_deco}
var
   i,j,k,l,mu,pp,o   : integer;
   ki,kj,li,lj	     : integer;
   b,bf		     : tsyml;
   bm		     : array[1..t_sym] of tsyml;
   {b,bf		     : array[1..t_sym] of integer;}
   b3		     : array[-1..1] of boolean;
   sg,o1,o0,oo,objnu : integer;
   reg,obreak	     : boolean;
   rgp		     : array[1..npth] of boolean;
   glist	     : array [1..t_sym] of noderec5;
   ngl		     : array [1..t_sym] of integer;
   nulist	     : array[1..t_sym] of i1list;
   vli,vlj	     : i5list;
   eps		     : myreal;
   iw5		     : i5;
   bfp,bmp	     : integer;
   function checkfli(x: i3):integer;
   Var
      i,answer : integer;
   begin
      answer:=0;
      if(nfli>0) then begin
	 i:=0; repeat i:=i+1;
	    if ((fli[i,1]=x[1]) and (fli[i,2]=x[2]) and (fli[i,3]=x[3])) then answer:=i;
	 until ((answer>0) or (i=nfli));
      end;
      if (answer=0) then begin
	 nfli:=nfli+1;
	 fli[nfli]:=x;
	 answer:=nfli;
      end;
      checkfli:=answer;
   end; { checkfli }
   procedure pth25d(np : integer;  x : vpath; Var y:i5list);
   Var i,o,pp : integer;
   begin
      for i:=1 to np do begin
	 for o:=1 to 5 do y[i,o]:=0;
	 if (i=1) then begin
	    pp:=x[1];
	    y[i]:=bs[pp];
	 end else begin
	    pp:=pp+x[i];
	    for o:=1 to 5 do y[i,o]:=y[i-1,o]+bs[po(pp),o];
	 end;
      end;
   end; { pth25d }
   function is_not_in_list(nn : integer;y:noderec5;z:i6):boolean;
   Var i,j,o   : integer;
      found,eq : boolean;
      d	       : i5;
   begin
      j:=0;
      is_not_in_list:=true; found:=false;
      repeat j:=j+1;
	 for o:=1 to 5 do d[o]:=y[j,o]-z[o]; {IGNORE the "height"!}
	 {Comparison: difference 5D vector d must be (i,i,i,i,i)=0 in real space}
	 eq:=true; for o:=2 to 5 do if (d[o]<>d[1]) then eq:=false;
	 if eq then found:=true;
      until ((j=nn) or found);
      if found then is_not_in_list:=false;
   end;
   procedure pathlist5(np : integer; x,xs: vpath; Var nn:integer;Var y:noderec5;io:integer);
   {x is path; xs is status; nn index of y, y is pathlist}
   Var
      i,j,o,pp : integer;
      y0       : i5;
      z0       : i6;
   begin
      for i:=1 to np do begin
	 if (i=1) then begin
	    pp:=x[1];
	    y0:=bs[pp];
	 end else begin
	    pp:=pp+x[i];
	    for o:=1 to 5 do y0[o]:=y0[o]+bs[po(pp),o];
	 end;
	 for o:=1 to 5 do z0[o]:=y0[o]; z0[6]:=io;
	 if (((xs[i]=1) or (xs[i]<0)) and is_not_in_list(nn,y,z0)) then begin {when status=1, AND the vertex not already in the list, add to the list}
	    nn:=nn+1;
	    for o:=1 to 6 do y[nn,o]:=z0[o];
	    if(xs[i]=1) then y[nn,7]:=0 else y[nn,7]:=xs[i];
	 end;
      end;
   end;
   procedure patlist5(x : i10;Var nn:integer; Var y:noderec5;io:integer);
   {Obtain list of 5D vectors from the "pattern record".
    This contains 10 integers indicating status of vertex along that
    direction; 0=empty, 1=occupied, 2=wildcard.}
   Var
      i,j,o : integer;
   begin
      for j:=1 to 10 do if ((x[j]=1) or (x[j]<0)) then begin
	 nn:=nn+1;
	 for i:=1 to 5 do y[nn,i]:=bs[j,i];
	 y[nn,6]:=io;
	 if(x[j]=1) then y[nn,7]:=0 else y[nn,7]:=x[j];
      end;
   end; { i10_2_v5 }
   function register5(nn :integer; a,b:noderec5;nua,nub:i1list;ca,cb:v5):boolean;
   Var
      i,j,k,o,io : integer;
      u,v	 : v5;
      w		 : v2;
      sum	 : myreal;
      ax,bx	 : array[1..listsize] of integer;
      regi,nofl	 : boolean;
   begin
      {writeln;}
      nofl:=true;
      for i:=1 to nn do begin ax[i]:=0; bx[i]:=0; end;
      for i:=1 to nn do begin
	 for j:=1 to nn do if ((bx[j]=0) and (ax[i]=0)) then begin
	    {if ((nua[i]=nub[j]) and (a[i,7]=b[j,7])) then begin}
	    if (a[i,7]=b[j,7]) then begin
	       for o:=1 to 2 do w[o]:=0.0;
	       for o:=1 to 5 do begin
		  u[o]:=a[i,o]-ca[o];
		  v[o]:=b[j,o]-cb[o];
		  for io:=1 to 2 do w[io]:=w[io]+(u[o]-v[o])*e[o,io];
	       end;
	       sum:=0.0;
	       for o:=1 to 2 do sum:=sum+w[o]*w[o];
	       {write(i:1,'-',j:1,' ');
		for o:=1 to 5 do write(u[o]:1:3,' ');
		write(' - ');
		for o:=1 to 5 do write(v[o]:1:3,' ');}
	       if (sum<eps) then begin
		  {write(' ',i:1,'+',j:1);}
		  ax[i]:=j; bx[j]:=i;
	       end;
	       {writeln('  ',ax[i]:1,' ',bx[j]:1);}
	    end;
	 end;
      end;
      regi:=true;
      k:=0; for i:=1 to nn do if (ax[i]=0) then regi:=false;
      {writeln(' ',regi);}
      register5:=regi;
   end; { register5 }
begin
   {writeln('TEST obj_name : ',obj_name[no]);}
   for i:=1 to n_sym do begin
      for o:=1 to 3 do obj_syml[no,i,o]:=1;
   end;
   eps:=1e-5; {for registering symmetry-related patterns}
   {First, generate all symmetry-related object-path constructs}
   objnu:=obj_nu[no,1];
   for i:=1 to n_sym do begin b[i]:=0;
      for j:=1 to t_sym do bm[i,j]:=0;
      obj_nu[no,i]:=objnu*defdp[i];
      obj_dpz[no,i]:=defdp[i];
      {writeln(i:2,obj_nu[no,i]:4,defdp[i]:4);}
      for o:=1 to 3 do obj_fla[no,i,o]:=obj_flad[no,o]*syml[i,o];
      {TEST THIS:}
      {for o:=1 to 3 do obj_fla[no,i,o]:=obj_flad[no,o];}
      {write('flad ',no:3,i:3,o:3,'  ');
      for o:=1 to 3 do write(obj_fla[no,i,o]:2); writeln;}
      sg:=1; {for path definition!!}
      if (syml[i,2]>1) then sg:=-1;
      {THIS vertex: pattern -> 5D}
      ngl[i]:=0;
      for o:=-io3d to io3d do if ((obj_node[no,o]=1) or (obj_node[no,o]<0)) then begin
	 ngl[i]:=ngl[i]+1;
	 for j:=1 to 5 do glist[i,ngl[i],j]:=0;
	 glist[i,ngl[i],6]:=o;
	 glist[i,ngl[i],7]:=obj_node[no,o];
      end;
      for j:=1 to 10 do begin
	 for o:=-io3d to io3d do begin
	    obj_pat_list[no,i,o,j]:=obj_pat[no,o,sym[i,j]];
	 end;
	 for o:=1 to obj_nipat[no] do begin
	    obj_ipat_list[no,i,o,j]:=obj_ipat[no,o,sym[i,j]];
	 end;
	 {set default 3D pattern anyway; 3D may be switched later on}
	 if (io3d=0) then for o:=-1 to 1 do if (o<>0) then
	    obj_pat_list[no,i,o,j]:=obj_pat[no,0,sym[i,j]];
	 for o:=1 to 2 do if (obj_bas[o,sym[i,j]]>0) then obj_bas_list[i,o]:=j;
      end;
      {write(obj_name[no],' ');
      for o:=-io3d to io3d do begin
      for j:=1 to 10 do write(obj_pat_list[no,i,o,j]:1); write(' ');
      end; writeln;}
      {PATTERN -> 5D}
      for o:=-io3d to io3d do PATLIST5(obj_pat_list[no,i,o],ngl[i],glist[i],o);
      if (obj_pth[no]>0) then begin
	 for j:=1 to 10 do if (sym[i,j]=obj_pth_def[no,1,1]) then
	    obj_pth_def[no,i,1]:=j;
	 for j:=2 to obj_pth[no] do
	    obj_pth_def[no,i,j]:=sg*obj_pth_def[no,1,j];
	 {for j:=1 to obj_pth[no] do write(obj_pth_def[no,i,j]:1,' ');
	 writeln('  obj_pth_def ',no:1,' ',i:1);}
      end;
      {PATH -> 5D}
      for o:=-io3d to io3d do PATHLIST5(obj_pth[no],obj_pth_def[no,i],obj_pth_stat[no,o],ngl[i],glist[i],o);
   end;
   for i:=1 to n_sym do begin
      for o:=1 to 5 do gor[i,o]:=0.0;
      for j:=1 to ngl[i] do begin
	 for o:=1 to 5 do iw5[o]:=glist[i,j,o];
	 nulist[i,j]:=dperp_offset(obj_nu[no,i],nui(iw5));
	 {dbg temp XXX dpmodel2}
	 {nulist[i,j]:=nui(iw5);}
	 if (center=1) then for o:=1 to 5 do gor[i,o]:=gor[i,o]+glist[i,j,o];
      end;
      if (center=1) then begin
	 for o:=1 to 5 do gor[i,o]:=gor[i,o]/ngl[i];

	 {for j:=1 to ngl[i] do begin
	 write('dbgXXX dpmod=',dpmodel:1,' ',defdp[i]:1,' ',obj_name[no],' ',i:2,' ',j:2,'   ');
	    for o:=1 to 5 do write(glist[i,j,o]:1,' ');
	 write('   ',glist[i,j,6]:1,' ',glist[i,j,7]:1);
	    write('   ',nulist[i,j]:1,'  ');
	    for o:=1 to 5 do write((glist[i,j,o]-gor[i,o]):1:2,' ');
	    writeln;
	 end;
	 write('	OBJ_ORIG : ');
	 for o:=1 to 5 do write(gor[i,o]:1:2,' ');
	 writeln;}
	 
      end;
   end;
   nfli:=0;
   for i:=1 to n_sym do begin
      for o:=1 to 3 do ofla[o]:=obj_fla[no,i,o];
      bf[i]:=checkfli(ofla);
      {writeln('dbg-fli ',obj_name[no],' ',i:3,bf[i]:5);}
   end;
   if(center=1) then begin
      {writeln(obj_name[no],' FLIDBG ',nfli:1);
      for i:=1 to nfli do begin
      write('fli ',i:1,'  ');
      for o:=1 to 3 do write(fli[i,o]:1,' '); writeln;
      end;}
   end;
   {write(obj_name[no],': ');}
   {Following is symmetry reduction of the object}
   for i:=1 to n_sym-1 do begin
      for j:=i+1 to n_sym do if (b[i]<=0) and (b[j]<=0) then begin
	 reg:=(ngl[i]=ngl[j]); {status+3D may cause inequality here}
	 if reg then reg:=register5(ngl[i],glist[i],glist[j],nulist[i],nulist[j],gor[i],gor[j]);
	 {finally, test for user-supplied object symmetry-breaking}
	 if (reg) then begin
	    obreak:=false;
	    for o:=1 to 3 do begin
	       if not(obj_fla[no,i,o]=obj_fla[no,j,o]) then begin
		  obreak:=true;
		  if not(obj_break[no]) then begin
		     writeln('OBJ-SYMM-BREAK: ',no:1,' ',obj_name[no]);
		     obj_break[no]:=true;
		  end;
	       end;
	    end;
	    {for o:=1 to 3 do write(obj_fla[no,i,o]:1,'\',obj_fla[no,j,o]:1,' '); writeln('   TEST OBJ-S-B ',obj_name[no]);}
	 end;
	 if (reg) then begin
	    {nfold is 5 or 10}
	    ki:=(i-1) div nfold; kj:=(j-1) div nfold;
	    {compute site symmetries}
	    if (ki=kj) then begin {5-fold symmmetry}
	       obj_syml[no,i,1]:=obj_syml[no,i,1]+1;
	    end;
	    if ((syml[i,3]=syml[j,3])) then begin {mirror}
	       obj_syml[no,i,2]:=obj_syml[no,i,2]+1;
	    end;
	    if (syml[i,2]=syml[j,2]) then begin {2-fold-s}
	       obj_syml[no,i,3]:=obj_syml[no,i,3]+1;
	    end;
	    b[j]:=i;
	 end;
      end; {j-loop through n_sym}
   end; {i-loop through n_sym}
   mu:=0; for i:=1 to n_sym do begin
      if (b[i]=0) then begin
	 bm[i,bf[i]]:=bm[i,bf[i]]+1;
	 bmp:=bm[i,bf[i]];
	 obj_bas_dp[i,bf[i],bmp]:=1;
	 obl[i,bf[i],bmp]:=obj_bas_list[i];
	 mu:=mu+1;
	 {writeln('dbg-obl ',i:3,bf[i]:3,bmp:4,obl[i,bf[i],bmp,1]:5,obl[i,bf[i],bmp,2]:3);}
      end else begin
	 pp:=b[i];
	 bfp:=bf[i];
	 bm[pp,bfp]:=bm[pp,bfp]+1;
	 bmp:=bm[pp,bfp];
	 {this is for grey objects like HBS bond}
	 obj_bas_dp[pp,bfp,bmp]:=defdp[b[i]]*defdp[i];
	 obl[pp,bfp,bmp]:=obj_bas_list[i];
	 {writeln('dbg-obl-2 ',pp:3,bfp:3,bmp:4,obl[pp,bfp,bmp,1]:5,obl[pp,bfp,bmp,2]:3);}
      end;
   end;
   mu:=0;
   for i:=1 to n_sym do if (b[i]=0) then begin
      mu:=mu+1;
      for j:=1 to nfli do
	 for k:=1 to bm[i,j] do begin
	    obl[mu,j,k]:=obl[i,j,k];
	    obj_bas_dp[mu,j,k]:=obj_bas_dp[i,j,k];
	    {writeln('obl-dbg-3 ',mu:3,j:3,k:3,obl[mu,j,k,1]:5,obl[mu,j,k,2]:3);}
	 end;
      {this must be 1, not  io3d! (Because 3d may be switched on later)}
      for o:=-1 to 1 do
	 obj_pat_list[no,mu,o]:=obj_pat_list[no,i,o];
      for o:=1 to obj_nipat[no] do begin
	 obj_ipat_list[no,mu,o]:=obj_ipat_list[no,i,o];
      end;
      obj_dpz[no,mu]:=obj_dpz[no,i];
      obj_nu[no,mu]:=obj_nu[no,i];
      obj_fla[no,mu]:=obj_fla[no,i];
      obj_syml[no,mu,1]:=obj_syml[no,i,1];
      for o:=2 to 3 do obj_syml[no,mu,o]:=obj_syml[no,i,o] div obj_syml[no,i,1];
      {pattern origin only for screw2}
      if (obj_pth[no]>0) then begin
	 for k:=1 to obj_pth[no] do
	    obj_pth_def[no,mu,k]:=obj_pth_def[no,i,k];
      end;
      if(center=1) then begin
	 gor[mu]:=gor[i];
	 {
	 write('DBG(i,nfli,bm,obl) ',obj_name[no],' ',i:1,' ',nfli:1,' ');
	 for j:=1 to nfli do write(bm[i,j]:1,' ');
	 write(' ');
	 for j:=1 to nfli do
	 for k:=1 to bm[i,j] do write(obl[i,j,k,1]:1,':',obl[i,j,k,2]:1,' ');
	 writeln;
	 }
      end;
   end;
   ob:=b;
   {for i:=1 to mu do begin
    for j:=1 to bm[i] do begin
    obj_bas_dp[i,j]:=obj_bas_dp[i,j]*obj_bas_dp[i,1];
      end;
   end;}
   {if (center=1) then for i:=1 to n_sym do if (b[i]=0 and bm[i]*mu<>n_sym) then begin
   writeln('Basis_multiplicity error: ',i:1,' ',bm[i]:1,' ',mu:1,' ',n_sym:1);
   end;}
   obj_mu[no]:=mu;
   obj_bas_mu:=bm[1,1]; {must be same for all flavours all orientations}
end; { union }
procedure obj_union(no : integer);
Var
   center,omu,i,j,o : integer;
begin
   {we call obj_union first without centering the object, then
   with centering. If centering increased symmetry, the centering
   shift "obj_ori" is exported.}
   for i:=1 to n_sym do for o:=1 to 2 do objori[no,i,o]:=0.0;
   center:=0;
   obj_union2(no,center,or5);
   omu:=obj_bas_mu;
   center:=1;
   obj_union2(no,center,or5);
   if (autocc[no]=1) then center:=1 else center:=0;
   if((omu<obj_bas_mu) and (center=0)) then begin
      writeln('mkdeco : autocentering would increase symmetry of ',obj_name[no],". Use 'autocenter' wihin object if needed");
   end;
   {writeln(obj_name[no],' ',omu:1,' ',obj_bas_mu:1);}
   {if centering increased symmetry, use centering shift}
   bobor[no]:=0;
   {writeln(obj_name[no],' ',omu:1,' ',obj_bas_mu:1);}
   {NOTE: condition to use obj_break for preserving centering is only provisional solution: actually we want the first condition, but then avoid force-sym-break feature to disable centering}
   {if ((omu<obj_bas_mu) or obj_break[no]) then begin}
   if (center=1) then begin
      bobor[no]:=1;
      {write('USING SHIFTED ORIGIN');}
      for i:=1 to obj_mu[no] do begin
	 {previosly used condition: "if (obj_syml[no,mu,3]>1) then"...}
	 r5r2(obj_bas_list[i],or5[i],objori[no,i]);
	 {write('DBG ',obj_name[no],' ',i:1,' ');
	 for o:=1 to 5 do write(or5[i,o]:1:5,' '); writeln;}
	 for o:=1 to 2 do or2[no,i,o]:=0.0;
	 for j:=1 to 5 do
	    for o:=1 to 2 do or2[no,i,o]:=or2[no,i,o]+or5[i,j]*e[j,o];
      end;
   end else begin
      for i:=1 to obj_mu[no] do begin
	 for o:=1 to 5 do or5[i,o]:=0;
	 for o:=1 to 2 do or2[no,i,o]:=0;
      end;
   end;
end;
procedure lift(j1,j2: integer;ox,x0: v2;Var dx:v2;b2:i2;bdp:integer;Var y5:v5);
{v are coefficients along k1 and k2 basis directions.
We lift real coeffs v and find integer representation
using inflation matrix infl[10,5] defined in tools.p}
Var
   a,b,nm,nm2,o,i,j,nu : integer;
   xpa,xpe,eps,cc      : myreal;
   amax,bmax	       : integer;
   done		       : array[1..2] of boolean;
   x5		       : i5;
   x		       : v2;
begin
   for o:=1 to 2 do x[o]:=x0[o]+bdp*ox[o];
   amax:=round(liftperplim)+1;
   bmax:=round(liftperplim/tau)+1;
   for o:=1 to 5 do y5[o]:=0.0;
   for o:=1 to 2 do done[o]:=false;
   nm:=0; nm2:=0;
   for o:=1 to 2 do if not(done[o]) then begin
      for a:=-amax to amax do
	 for b:=-bmax to bmax do begin
	    xpa:=a*tau+b;
	    xpe:=a-b*tau;
	    if (xpe<liftperplim) then begin
	       if (abs(xpa-x[o])<lifteps) then begin
		  dx[o]:=x[o]-xpa;
		  nm:=nm+1;
		  for i:=1 to 5 do
		     y5[i]:=y5[i]+a*infl[b2[o],i]+b*bs[b2[o],i];
		  done[o]:=true;
	       end;
	    end;
	 end;
   end;
   if (nm=2) then begin
      for o:=1 to 5 do x5[o]:=round(y5[o]);
      nu:=nui(x5);
      for o:=1 to 5 do y5[o]:=x5[o];
      {write('R_atom = (',v[1]:1:5,', ',v[2]:1:5,')');
      write('  5D: '); for o:=1 to 5 do write(y5[o]:3); writeln;}
   end else begin
      {Try half-integer lift}
      amax:=2*round(liftperplim)+1;
      bmax:=2*round(liftperplim/tau)+1;
      for o:=1 to 5 do y5[o]:=0.0;
      for o:=1 to 2 do done[o]:=false;
      nm2:=0;
      for o:=1 to 2 do if not(done[o]) then begin
	 for a:=-amax to amax do
	    for b:=-bmax to bmax do begin
	       xpa:=(a*tau+b)/2.0;
	       xpe:=(a-b*tau)/2.0;
	       if (xpe<liftperplim) then begin
		  if (abs(xpa-x[o])<lifteps) then begin
		     dx[o]:=x[o]-xpa;
		     nm2:=nm2+1;
		     {for i:=1 to 5 do write(y5[i]:1:1,' '); writeln;}
		     for i:=1 to 5 do
			y5[i]:=y5[i]+(a*infl[b2[o],i]+b*bs[b2[o],i])/2.0;
		     done[o]:=true;
		  end;
	       end;
	    end;
      end;
      if ((nm<2) and (nm2<2)) then begin
	 write(obj_name[j1],', ORBIT ',j2:1,' -> ');
	 writeln('Couldnt lift R_atom = (',v[1]:1:5,', ',v[2]:1:5,').');
	 writeln('Increase perpliftlim, or increase lifteps, or disable lifting.');
	 halt;
      end;
      if ((nm>2) or (nm2>2)) then begin
	 write(obj_name[j1],', ORBIT ',j2:1,' -> ');
	 writeln('Lift is ambiguous, decrease lifteps or perpliftlim.');
	 halt;
      end;
      if ((nm=2) and (nm2=2)) then begin
	 write(obj_name[j1],', ORBIT ',j2:1,' -> ');
	 writeln('Ambiguous lift: match for integer and half-integer lift.');
	 halt;
      end;
      if ((nm2=2) and (first_atom)) then begin
	 write('Object ',obj_name[j1],', orbit ',j2:1,': ');
	 {for o:=1 to 5 do write(y5[o]:1:1,' ');}
	 writeln('half-integer lift (OK).');
      end;
   end;
end; { lift }

procedure READ_RAW_DECO;
var
   i,j,k,l,m,n,o    : integer;
   txt,txo,txdbg    : text;
   ch		    : char;
   ro		    : v2list;
   ro5		    : v5list;
   bo		    : atomreclist;
   sline	    : mystring;
   sll		    : string[250];
   kl,jo,ji,im,jato : integer;
   b2		    : i2;

   function obj_label(x :integer):integer; {in read_deco_human}
   var i,j,o : integer;
      b	     : boolean;
   begin     
      i:=0; repeat i:=i+1; b:=(obj_code[i]=x) until ((i=n_obj) or b);
      if b then obj_label:=i else obj_label:=0;
   end; { obj_label }
begin
   rewrite(txdbg,'debug.dat');
   blift:=false;
   io3d:=0;
   if (figout) then writeln(' -> INPUT/objects.fig') else writeln(' -> INPUT/objects.def');
   if (figout) then rewrite(txo,'INPUT/objects.fig') else rewrite(txo,'INPUT/objects.def');
   reset(txt,'INPUT/symmetry.def');
   readln(txt,n_sym,nfold,nmir,n2s);
   for i:=1 to n_sym do begin
      for j:=1 to 10 do read(txt,sym[i,j]);
      read(txt,defdp[i]);
      for o:=1 to 3 do read(txt,syml[i,o]);
      readln(txt);
   end;
   close(txt);
   for i:=1 to t_obj do
      for j:=1 to t_sym do obj_at_mu[i,j]:=0;
   for i:=1 to t_orb do begin mu_sb[i]:=0; end;
   ncode:=0; n_obj:=0; n_orb:=0; nstar:=0;
   n_pxyz:=0; dpmodel:=1;
   for i:=1 to t_obj do
      for o:=1 to 3 do obj_flad[i,o]:=0; {reset object symbreak status}
   for i:=1 to t_obj do obj_break[i]:=false;
   reset(txt,'INPUT/mkdeco.inp');
   readln(txt,deconame);
   writeln(txo,deconame);
   readln(txt,o); figout:=(o=1);
   readln(txt,nu_offs); {nu offset}
   writeln(txo,nu_offs:2,'  # NU_OFFSET (discrete perp-space)');
   writeln('Recentering discrete perp: ',nu_offs:1);
   readln(txt,dpmodel);
   write(txo,dpmodel:2,'  # dpmodel: ');
   if (dpmodel=0) then writeln('DPmodel==0: flavours OFF');
   if (dpmodel=1) then writeln('DPmodel==1: flavours ON');
   if (dpmodel=2) then writeln('DPmodel==2: flavours even-odd');
   if (dpmodel=0) then writeln(txo,' flavours OFF (0)');
   if (dpmodel=1) then writeln(txo,' flavours ON (1)');
   if (dpmodel=2) then writeln(txo,' flavours even-odd (2)');
   defdpz[1,1]:=1; {dp-mirror}
   defdpz[1,2]:=0; {dp-screw}
   defdpz[0,1]:=1;
   defdpz[0,2]:=0;
   readln(txt,defdpz[-1,1],defdpz[-1,2]);
   read(txt,o); blift:=(o=1); repre5:=blift;
   if(blift)then begin readln(txt,liftperplim,lifteps);
      writeln('Lifting using perp_limit ',liftperplim:1:2,' and tolerance ',lifteps:1:5);      
   end else readln(txt);
   write(txo,o:2,' ');
   if(blift) then writeln(txo,liftperplim:1:2,' ',lifteps:1:5,'  # perp-lift: [1|0] [PERP-BOUND] [TOL]') else writeln(txo,' # no 5D lift');
   readln(txt,o); deco:=(o=1); {set definition type atoms/tile-ham}
   if deco then write('Input type objects-atoms from decofile ');
   if not(deco) then write('Input type tile-hamiltonian from ');
   writeln(txo,o:2,'  # decotype: 0=tile-ham; 1=atomic decoration; -1=undef');
   writeln(deconame);
   readln(txt,io3d);
   if(io3d=1)then writeln('Object definition with 3D-extension.');
   writeln(txo,io3d:2,'  # 3D extension of object definitions');
   writeln('-------------------------------------------------------');
   readln(txt,ntobj);
   writeln(txo,ntobj:2,'  # total defined objects');
   writeln(txo,'# ---------------- start objects defs ------------------');
   if (ntobj>t_obj) then begin
      writeln('Insufficient array dim t_obj=',t_obj:1,'. Decofile requires ',ntobj:1,'.');
   end;
   for n_obj:=1 to ntobj do begin
      readln(txt,obj_name[n_obj]);
      readln(txt,obj_node[n_obj,0]);
      readln(txt,nobcolors[n_obj]);
      {writeln(txdbg,n_obj:1,' ',nobcolors[n_obj]:1,' debug nobcolors 1');}
      readln(txt,autocc[n_obj]);
      readln(txt,obj_nu[n_obj,1]);
      if (dpmodel=0) then obj_nu[n_obj,1]:=0;
      {read(txt,obj_ninst[n_obj]);
      if(obj_ninst[n_obj]>0) then for o:=1 to obj_ninst[n_obj] do read(txt,obj_inst[n_obj,o]);
      readln(txt);}
      for o:=1 to 3 do read(txt,obj_flad[n_obj,o]); readln(txt);
      for i:=1 to 2 do
	 for j:=1 to 10 do obj_bas[i,j]:=0;
      read(txt,o); obj_bas[1,o]:=1;
      read(txt,o); obj_bas[2,o]:=1;
      readln(txt);
      readln(txt,obj_nipat[n_obj]);
      k:=obj_nipat[n_obj];
      if(k>0) then begin
	 for i:=1 to k do begin
	    read(txt,obj_inst[n_obj,i]);
	    for o:=1 to 10 do read(txt,obj_ipat[n_obj,i,o]);
	 end;
	 readln(txt);
      end;
      read(txt,obj_nipth[n_obj]);
      for i:=1 to obj_nipth[n_obj] do read(txt,obj_ipth[n_obj,i]);
      readln(txt);
      for o:=1 to 10 do read(txt,obj_pat[n_obj,0,o]); readln(txt);
      for o:=-1 to 1 do if (o<>0) then
	 obj_pat[n_obj,o]:=obj_pat[n_obj,0];
      read(txt,obj_pth[n_obj]);
      for i:=1 to obj_pth[n_obj] do begin
	 read(txt,obj_pth_def[n_obj,1,i]);
	 for o:=-1 to 1 do obj_pth_stat[n_obj,o,i]:=1;
      end;
      for i:=1 to obj_pth[n_obj] do begin
	 read(txt,obj_pth_stat[n_obj,0,i]);
      end;
      readln(txt); {end path definition}
      if(io3d>0)then begin
	 readln(txt,obj_node[n_obj,-1],obj_node[n_obj,1]);
	 for o:=-io3d to io3d do if (o<>0) then begin
	    for i:=1 to 10 do read(txt,obj_pat[n_obj,o,i]);
	 end;
	 readln(txt);
	 if(obj_pth[n_obj]>0) then begin
	    for o:=-io3d to io3d do if (o<>0) then begin
	       for i:=1 to obj_pth[n_obj] do read(txt,obj_pth_stat[n_obj,o,i]);
	    end;
	    readln(txt);
	 end;
      end;
      obj_code[n_obj]:=999;
      OBJ_UNION(n_obj); {inits the object}
      if (deco) then begin
	 jat:=0;
	 readln(txt,ntorb);
	 writeln('object ',obj_name[n_obj],', ',ntorb:1,' orbits ',(n_orb+1):1,'-',(n_orb+ntorb):1);
	 for jo:=1 to ntorb do begin
	    n_orb:=n_orb+1;
	    read(txt,x2[1],x2[2],v0[1],v0[2],zz,jch);
	    read(txt,ocolor); {user-defined object-color as an attribute of orbit}
	    read(txt,mu_sb[n_orb]);
	    if(mu_sb[n_orb]>0) then for o:=1 to mu_sb[n_orb] do read(txt,ch_sb[n_orb,o]);
	    readln(txt);
	    if (x2[1]=x2[2]) then begin
	       writeln('Wrong orbit def (i1==i2).'); halt;
	    end;
	    v:=v0;
	    for o:=1 to 10 do x10[o]:=0.0;
	    for o:=1 to 2 do x10[x2[o]]:=v[o];
	    k:=n_obj;
	    v10v2(obj_bas[1],obj_bas[2],x10,v);
	    bv1:=(abs(v[1])>0.0001); 
	    bv2:=(abs(v[2])>0.0001);
	    bv12:=(abs(v[2]-v[1])<0.0001);
	    i:=n_pxyz; {store}
	    {writeln('TTTA ',v[1]:1:3,' ',v[2]:1:3,' ',objori[n_obj,1,1]:1:3,' ',objori[n_obj,1,2]:1:3);}
	    for o:=1 to 2 do v[o]:=v[o]-objori[n_obj,1,o];
	    for o:=1 to 2 do begin
	       pxyz2orb[n_pxyz+o]:=n_orb;
	       pxyz[n_pxyz+o]:=v[o];
	       pxyz_rsd[n_pxyz+o]:=0.0;
	       papo[o]:=n_pxyz+o;
	    end;
	    n_pxyz:=n_pxyz+2;

	    {mirror treatment}
	    if ((obj_syml[n_obj,1,2]>1) and bv12) then begin
	       n_pxyz:=n_pxyz-1;
	       papo[2]:=papo[1];
	    end;
	    
	    {rotational symmetry constraint: when v=(0,0), no degrees
	    of freedom in deca plane.}
	    if ((obj_syml[n_obj,1,1]>1) and (absv(v)<0.0001)) then begin
	       n_pxyz:=i; {restore old value}
	       papo[1]:=0; papo[2]:=0;
	    end;
	    {here should go z-mirror, free unless z-mirror... (unimplemented)}
	    n_pxyz:=n_pxyz+1;
	    pxyz2orb[n_pxyz]:=n_orb;
	    papo[3]:=n_pxyz;
	    pxyz[n_pxyz]:=zz;
	    pxyz_rsd[n_pxyz]:=0.0;

	    { ZMIRROR stuff continued:
	     if(zmirror=1) then begin
	     zmi:=bound(2*zmirror-zz);
	     if(abs(zmi-zz)<0.0001)then begin
	     papo[3]:=0;
	       end else begin
	     n_pxyz:=n_pxyz+1;
	     pxyz2orb[n_pxyz]:=n_orb;
	     papo[3]:=n_pxyz;
	     pxyz[n_pxyz]:=zz;
	     pxyz_rsd[n_pxyz]:=0.0;
	       end;
	    end;
	     }
	    
	    {write('XXX ',k:1,' ',obj_bas_mu:1,' ',obj_mu[k]:1,' ');
	    writeln(obj_name[n_obj]);}
	    for i:=1 to obj_mu[n_obj] do begin
	       jat:=0; jato:=0;
	       for ji:=1 to nfli do begin
		  for l:=1 to obj_bas_mu do begin
		     jat:=jat+1;
		     first_atom:=((l=1) and (i=1));
		     if blift then lift(n_obj,n_orb,objori[n_obj,i],v,dv,obl[i,ji,l],obj_bas_dp[i,ji,l],ro5[jat]);
		     if ((papo[1]>0) and (papo[2]>0)) then begin
			if (papo[1]<>papo[2]) then begin
			   pxyz_rsd[n_pxyz-2]:=dv[1];
			   pxyz_rsd[n_pxyz-1]:=dv[2];
			end else pxyz_rsd[n_pxyz-1]:=dv[1];
		     end else begin
			if (papo[1]>0) then pxyz_rsd[n_pxyz-1]:=dv[1];
			if (papo[2]>0) then pxyz_rsd[n_pxyz-1]:=dv[2];
		     end;
		     for o:=1 to 3 do bo[jat,o]:=papo[o];
		     for n:=1 to 2 do
			bo[jat,3+n]:=obl[i,ji,l,n];
		     bo[jat,6]:=defdpz[obj_dpz[n_obj,i],1];  {dp-mirror}
		     bo[jat,7]:=defdpz[obj_dpz[n_obj,i],2];  {dp-screw}
		     bo[jat,8]:=0; {reserved for z-mirror}
		     bo[jat,9]:=(ocolor-1)*nfli+ji; {flavour inherited from object}
		     {THIS IS SPECIAL CASE OF GREY OBJECT LIKE HBS_BOND}
		     if (obj_bas_dp[i,ji,l]<0) then begin
			bo[jat,6]:=defdpz[-obj_dpz[n_obj,i],1];
			bo[jat,7]:=defdpz[-obj_dpz[n_obj,i],2];
		     end;
		     if (bo[jat,7]<0) then bo[jat,7]:=0;
		     for m:=1 to 2 do ro[jat,m]:=0.0;
		     for n:=1 to 2 do b2[n]:=obl[i,ji,l,n];
		     for m:=1 to 2 do for n:=1 to 2 do
			ro[jat,m]:=ro[jat,m]+v[n]*re[b2[n],m];
		  end; {end l-loop through identical (registered) object sets}
		  {write(jato:1,' ',jat:1);}
		  ORB_UNION(jato+1,jat,ro,ro5,bo);
		  {writeln(' x ',jato:1,' ',jat:1);}
		  {writeln('ORBUNION ',jato:1,' ',jat:1,' ',iat:1,' ',ji:1,' ');}
		  if ((mu_sb[n_orb]>0) and (jat-jato<>mu_sb[n_orb])) then begin
		     writeln('CH-Sym-Break: obj orb ',n_obj:1,' ',n_orb:1,' sb_mu=',mu_sb[n_orb]:1,' incompatible with the calculated mu=',(jat-jato):1,'.');
		     {halt;}
		  end;
		  jato:=jat;
	       end; {end ji-loop through object flavours}
	       {writeln(k:5,i:5,obj_bas_mu:5,iat:5);}
	       m:=n_orb;
	       for kl:=1 to jat do begin
		  obj_at_mu[k,i]:=obj_at_mu[k,i]+1;
		  n:=obj_at_mu[k,i];
		  ra[k,i,n]:=ro[kl];
		  rz[k,i,n]:=zz+obj_dpz[k,i]*0.5;
		  ra5[k,i,n]:=ro5[kl];
		  bi[k,i,n]:=bo[kl];
		  {test}
		  if blift then begin
		     for o:=1 to 2 do v0[o]:=0.0;
		     for o:=1 to 2 do
			for j:=1 to 5 do v0[o]:=v0[o]+ro5[kl,j]*e[j,o];
		     cmp:=dist(ro[kl],v0);
		     if (cmp>0.0001) then begin
			writeln('Lift error: 	');
			write('REAL: ');
			for o:=1 to 2 do write(ro[kl,o]:1:5,' ');
			write('  LIFTED: ');
			for o:=1 to 2 do write(v0[o]:1:5,' ');
			writeln;
			halt;
		     end;
		  end;
		  if (mu_sb[n_orb]>0) then begin
		     atp[k,i,n,1]:=ch_sb[n_orb,kl];
		  end else atp[k,i,n,1]:=jch; {chemistry}
		  atp[k,i,n,2]:=n_orb; {absolute orbit #}
		  atp[k,i,n,3]:=jat;
	       end; {end kl-loop through individual atoms}
	    end; {end i-loop through object multiplicity}
	 end; {end jo-loop through orbits bound to an object}
      end; {end if-deco}
      if not(deco) then begin
	 read(txt,obj_ham[n_obj],o);
	 if(o>0)then begin {preset a desired object density - former "+" prefix}
	    ncode:=ncode+1;
	    read(txt,obj_stx[n_obj]);
	    obj_code[n_obj]:=ncode;
	 end;
	 readln(txt);
	 {writeln(obj_name[n_obj],' ',ncode:1,' ',obj_code[n_obj]:1);}
      end;
      readln(txt,sline); {the extra line}
      writeln(txo,n_obj:1,' ',obj_code[n_obj]:1,' ',obj_name[n_obj],' # object code&name');
      if (deco) then writeln(txo,jat:2,'  # atoms bound to object') else begin
	 if (obj_code[n_obj]<99) then writeln(txo,obj_stx[n_obj]:1:6,' ',obj_ham[n_obj]:1:5,' # obj_chp obj_ham') else writeln(txo,obj_ham[n_obj]:1:5,' # obj_ham');
      end;
      {write(txo,obj_ninst[n_obj]:2,'  ');
      if(obj_ninst[n_obj]>0) then for o:=1 to obj_ninst[n_obj] do write(txo,' ',obj_inst[n_obj,o]:1);
      writeln(txo,' # flip_node_instances');}
      write(txo,obj_nipat[n_obj]:2,' ');
      if(obj_nipat[n_obj]>0) then for k:=1 to obj_nipat[n_obj] do write(txo,obj_inst[n_obj,k]:1,' ');
      writeln(txo,'  # flip_pattern instances');
      write(txo,obj_nipth[n_obj]:2,' ');
      if(obj_nipth[n_obj]>0) then for k:=1 to obj_nipth[n_obj] do write(txo,obj_ipth[n_obj,k]:1,' ');
      writeln(txo,'  # flip_boundary instances');
      
      writeln(txo,obj_mu[n_obj]:2,'  # object multiplicity');
      writeln(txo,obj_pth[n_obj]:2,'  # object "PATH" definition length');
      if(nfli>1) then begin
	 writeln(txo,nfli:2,'  # colors via symmetry-breaking');
	 {for k:=1 to nfli do begin
	    for o:=1 to 3 do write(txo,fli[k,o]:1,' ');
	    if(k>1) then writeln(txo) else writeln(txo,'  # ',nfli:1,' lines: obj-sym-break: ROT(5|10),MIR,ROT(2+screw|dpinv)');
	 end;}
      end else writeln(txo,' 1  # no symmetry breaking (single object flavour)');
      if(nfli>1) then begin
	 for k:=1 to obj_mu[n_obj] do begin
	    for l:=1 to nfli do
	       write(txo,obl[k,l,1,1]:1,' ',obl[k,l,1,2]:1,'  ');
	    if(k>1) then writeln(txo) else writeln(txo,' # ',obj_mu[n_obj]:1,' lines: b1 b2 ...');
	 end;
      end;
      {writeln(txdbg,n_obj:1,' ',nobcolors[n_obj]:1,' debug nobcolors 2');}
      writeln(txo,nobcolors[n_obj]:2,'  # user-defined object colors');
      if (bobor[n_obj]=0) then writeln(txo,' 0  # no object centering') else writeln(txo,' 1  # object centering');
      if(bobor[n_obj]>0) then begin
	 for k:=1 to obj_mu[n_obj] do begin
	    for o:=1 to 5 do write(txo,or5[k,o]:1:2,' '); if(k>1) then writeln(txo) else writeln(txo,'  # ',obj_mu[n_obj]:1,' lines: object origin');
	 end;
      end;
      writeln(txo,obj_node[n_obj,0]:2,' # root-node status or instance');
      if (io3d>0) then begin
	 for o:=-io3d to io3d do if (o<>0) then write(txo,obj_node[n_obj,o]:1,' ');
	 writeln(txo,' # 3d extension for "root node": is vertex below&above present?');
      end;
      {obj_pth_status is the same for all symmetries}
      for m:=1 to obj_mu[n_obj] do begin
	 write(txo,obj_nu[n_obj,m]:1,'  ');
	 for o:=-io3d to io3d do begin
	    for k:=1 to 10 do write(txo,obj_pat_list[n_obj,m,o,k]:1,' ');
	    write(txo,'  ');
	 end;
	 for o:=1 to obj_nipat[n_obj] do begin
	    for k:=1 to 10 do write(txo,obj_ipat_list[n_obj,m,o,k]:1,' ');
	    write(txo,'  ');
	 end;
	 {if (m=1) then begin write(txo,'  # ',obj_mu[n_obj]:1,' lines: pattern ');
	    for o:=-io3d to io3d do write(txo,o:1,' ');
	 end;}
	 writeln(txo);
      end;
      if (obj_pth[n_obj]>0) then begin
	 for o:=-io3d to io3d do begin
	    for m:=1 to obj_pth[n_obj] do write(txo,obj_pth_stat[n_obj,o,m]:1,' ');
	    write(txo,' ');
	 end;
	 write(txo,'  # 1 line: pth_stat ');
	 for o:=-io3d to io3d do write(txo,o:1,' ');
	 writeln(txo);
	 for m:=1 to obj_mu[n_obj] do begin
	    for k:=1 to obj_pth[n_obj] do 
	       write(txo,obj_pth_def[n_obj,m,k]:1,' '); write(txo,' ');
	    if (m=1) then writeln(txo,'  # ',obj_mu[n_obj]:1,' lines: pth_def') else writeln(txo);
	 end;
      end;
      writeln(txo,'# -------- end_of_object');
   end; {end n_obj-loop (reading objects)}
   close(txt);
   close(txo);
   n_obj:=ntobj;
end; { read_deco_human }
   
BEGIN
   BASIS;
   READ_RAW_DECO; {reads symmetry.def and decoration}
   for i:=1 to n_obj do begin
      ptr[i]:=i;
      ptra[i,1]:=1.0; ptra[i,2]:=1.0; ptrz[i]:=1.0;
      for o:=1 to 5 do ptra5[i,o]:=1.0;
   end;
   if (n_obj>0) then begin
   end;
   if (deco) then begin
      n_pxyz0:=n_pxyz;
      {determine orbit chemistries}
      for i:=1 to 3*t_orb do nn_pxyz[i]:=0;
      for i:=1 to n_obj do begin
	 for m:=1 to obj_mu[i] do begin
	    for n:=1 to obj_at_mu[i,m] do begin
	       add_chem_para(bi[i,m,n],atp[i,m,n,1]);
	    end;
	 end;
      end;
      {add "real2perp" shift to all dpxyzs}
      for i:=1 to n_pxyz do
	 for j:=1 to nn_pxyz[i] do dpxyz[i,j]:=pxyz_rsd[i];
      writeln(n_orb:1,' atomic orbit(s), ',n_pxyz:1,' positional parameter(s) -> INPUT/ap.def');
      rewrite(txt,'INPUT/ap.def');
      writeln(txt,n_pxyz:1,'  # positional parameters');
      for i:=1 to n_pxyz do begin
	 write(txt,i:2,' ',pxyz2orb[i]:2,' ',(pxyz[i]-pxyz_rsd[i]):1:14,'  ',nn_pxyz[i]:2,' ');
	 for j:=1 to nn_pxyz[i] do write(txt,ch_pxyz[i,j]:2,' ',dpxyz[i,j]:1:5,'  ');
	 writeln(txt);
      end;
      {object origins}
      n_obj_ori:=0;
      for i:=1 to n_obj do if (bobor[i]>0) then
	 for m:=1 to obj_mu[i] do n_obj_ori:=n_obj_ori+1;
      writeln(txt,n_obj_ori:1,'  # tiling_node -> object_center shifts');
      j:=0;
      for i:=1 to n_obj do if (bobor[i]>0) then
	 for m:=1 to obj_mu[i] do begin
	    j:=j+1;
	    for o:=1 to 2 do obj_ori[j,o]:=or2[i,m,o];
	    writeln(txt,j:1,' ',i:1,'  ',m:1,'   ',obj_ori[j,1]:1:14,' ',obj_ori[j,2]:1:14,'  ',obj_name[i]);
	 end;
      n_obj_ori:=0;
      for i:=1 to n_obj do begin
	 for m:=1 to obj_mu[i] do begin
	    if (bobor[i]>0) then n_obj_ori:=n_obj_ori+1;
	    writeln(txt,i:1,' ',m:1,' ',obj_at_mu[i,m]:1,' ',obj_name[i]);
	    for n:=1 to obj_at_mu[i,m] do begin
	       for j:=1 to 3 do write(txt,bi[i,m,n,j]:2,' ');
	       write(txt,'  ');
	       for j:=1 to 2 do write(txt,bi[i,m,n,j+3]:2,' ');
	       write(txt,'  ');
	       for j:=1 to 2 do write(txt,bi[i,m,n,j+5]:2,' ');
	       write(txt,' ');
	       if (bobor[i]>0) then write(txt,n_obj_ori:2,' ') else write(txt,' 0 ');
               write(txt,bi[i,m,n,9]:2);
	       {write(txt,"  1");} {standing for corrected ...}
	       write(txt,'   ');
	       for o:=1 to 2 do write(txt,atp[i,m,n,o]:2,' ');
	       write(txt,'   ');
	       if blift then for o:=1 to 5 do write(txt,ra5[i,m,n,o]:5:1,' ');
	       writeln(txt);
	    end;
	 end;
      end;
      close(txt);
      rewrite(txt,'INPUT/olab.def');
      writeln(txt,n_orb:1,' ',dpmodel:1,'  # of orbits; dpmodel; [orbit chem site_mu object]');
      for i:=1 to n_obj do
	 for j:=1 to obj_at_mu[i,1] do begin
	    if (j=1) or ((j>1) and (atp[i,1,j-1,2]<>atp[i,1,j,2])) then begin
	       o:=atp[i,1,j,2]; {orbit #}
	       write(txt,o:2,atp[i,1,j,1]:4,atp[i,1,j,3]:4,' ',i:1,' ');
	       write(txt,obj_name[i],' ');
	       if (mu_sb[o]>0) then begin
		  write(txt,' CH-Sym-Break: ');
		  for jk:=1 to mu_sb[o] do write(txt,ch_sb[o,jk]:1,' ');
	       end;
	       writeln(txt);
	    end;
	 end;
      close(txt);
   end;
END.


